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Abstract 

In recent years intercalated and pillared graphitic systems have come under increasing scrutiny 
because of their potential for modern energy technologies. While traditional ab initio methods such 
as the LDA give accurate geometries for graphite they are poorer at predicting physicial properties 
such as cohesive energies and elastic constants perpendicular to the layers because of the strong 
dependence on long-range dispersion forces. 'Stretching' the layers via pillars or intercalation 
further highlights these weaknesses. We use the ideas developed by [J. F. Dobson et al, Phys. Rev. 
Lett. 96, 073201 (2006)] as a starting point to show that the asymptotic C'^D^^ dependence of the 
cohesive energy on layer spacing D in bigraphene is universal to all graphitic systems with evenly 
spaced layers. At spacings appropriate to intercalates, this differs from and begins to dominate the 
C4^D~^ power law for dispersion that has been widely used previously. The corrected power law 
(and a calculated C3 coefficient) is then unsuccesfully employed in the semiempirical approach of 
[M. Hasegawa and K. Nishidate, Phys. Rev. B 70, 205431 (2004)] (HN). A modified, physicially 
motivated semiempirical method including some C^D"^ effects allows the HN method to be used 
successfully and gives an absolute increase of about 2 — 3% to the predicted cohesive energy, while 
still maintaining the correct C^D^^ asymptotics. 
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I. INTRODUCTION 



The graphite form of carbon is a discretely layered material. The sp^ hybridised orbitals 
keep the layers in a rigid hexagonal pattern while the tTj, orbitals help bind the layers. This 
weak interlayer binding gives graphite a small elastic constant (C33) perpendicular to the 
plane which allows graphite to be 'stretched' by pillaring (see eg. ref. ]]]) and intercalation 
(see eg. ref. [2|]) by other substances with potentially useful applications for Hydrogen storage 
and other new energy technology. 

Standard density functional theory (DFT)^ based approaches such as the LDA and GGA- 
are known (see ref. ^ for a summary) to have problems predicting the interlayer binding 
energy and interlayer elastic constant of graphite at its experimental layer separation. This 
is presumed to be caused by the inability of these functionals to accurately include the 
long-range London dispersion forces (often denoted van der Waals forces in DFT papers, a 
notation we adopt to maintain consistency with other work). LDA/GGA correspondingly 
predict an exponentially decreasing binding energy for D ^ Dq (where D is the interlayer 
separation distance and Dq = 3.337A is the experimental interlayer separation distance) as 
opposed to the correct power law behaviour. 

Various authors^'^'^i^'^iSiii^ have proposed corrections to the LDA/GGA results that yield 
an additional long-range attractive layer-layer potential of the form C^D^"^. By contrast 
Dobson, White and Rubio (DWR)^- have shown that the asymptotic power law behaviour 
for bigraphene is C^D'^ due to its unusual bandstructure near the K point,— suggesting 
that even these ab initio and semiempirical corrections to LDA/GGA miss some important 
physics. 

In this work we first show that the C^D^^ power law is universal to many-layered graphitic 
systems with uniform interlayer separation, including those with an infinite number of layers 
such as rare gas intercalated^, or pillared graphite. 

We then use our energy expression to calculate the correct C3 coefficient for bulk graphite 
and adapt the method of Hasegawa and Nishidate (HN)^ to emply a corrected power law, 
thereby permitting empirical modelling of the non-asymptotic region when D ^ Dq. This 
investigation suggests that the different power-law and coefficient could have effects on 
semiempirical and other methods which assume a C^D'^ decay of the dispersion potential 
but that such effect may dominate only in for D > Dq. 
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II. ASYMPTOTIC POWER LAW 



The success of the random-phase approximation (RPA) in generating a correlation en- 
ergy functional through the Adiabatic Connection Formula and Fluctuation-Dissipation 
Theorem (ACFFDT) with the correct power law for long-range dispersion forces is well 
studied.— For the case of graphene compounds, DWR— used a long- 
wavelength approximation to the bare density-density response (xo) function of graphene to 
prove a C^D~^ dispersion potential for bigraphene while also reproducing known results for 
other materials through the same method. 

If we assume (as in DWR) that the in-plane response of a graphene plane can be approx- 
imated for low surface-parallel wavenumber (g) by a homogenous system of similar physics 
then we can write the RPA equation for the interacting density-density response (x) as 
follows: 

X\{,q,z,z']u) = Xoiq,z,z';u) 

+ \ dxdyxo{q,z,x;u)w{q,x,y)xx{q,y,z';u) (1) 



where the integrals are one- dimensional and A is a coupling constant to be used in the 
adiabatic connection formula. In the case of a layered system where each layer is highly 
localised in z space and separated by a distance D so that Xo{q^ z, z'\ u) = J2f=io^ xiq'-, u)5{z — 
z')5{z — iD) we may rewrite equation ([1]) as a tensor equation over layer indices i and j 

Xx{q^U]D) =x{q,u)l 

+ Ax(g, u)w{q)n{qD)xx{q, u; D) (2) 



where w{q) = and [0]jj = cUi^j = e "^"^l* ■'I (0 < i,j < N) so that Xxiq^ z, z';u) = 
j:^^[x{q,u;D)],,5{z-tD)S{z'-jD). 

We can use the ACFFDT to write the correlation energy per layer of a two-dimensionally 
homogeneous system as 

^ /»1 /»oo /"OO poo poo 

Ec = :r dX du qdq dz dz' 

4vr^ Jo Jo Jo J^oo J~oo 

X [Xx{q, z, z'\ u) - x{q, z, z'; u)] w(g)e"^l^'^^l. (3) 

Remembering that dispersion comes entirely from inter-layer correlation effects and making 
use of the delta functions thus lets us calculate the energy per unit area per layer of an 



A^-layered system through 



h 



1 fOO fOO 



U^dw = - J J J 

X [J^xiq, u] D) - J^xiq, u; oo)] (4) 

where 

J'x{q,u-D) =w{q)j^TT[xx{q,u;D)n{qD)]. (5) 

Due to the high level of symmetry ft takes the form of a Toeplitz matrix. This allows us 
to make use of Szego's Theorem (ref. {23] contains a good review of Szego's Theorem and 
its applications) to calculate the trace in the limit ^ 00 (these equations can also be 
obtained by Fourier methods). Defining 

no = t ^.e-^ = (6) 

cosh(giJ) — cos(4) 

as the Fourier transform of the tensor elements of Q, we then find 

Tx{q,u-D) = ^i:i[{i-\cn)-'c^i\ 
1 f\^_ cm 



sinh{qD)C 



(7) 



v/[cosh(gD) - ACsinh(gL))]2 - 1 

where C = x{q, u)w{q). 

For stretched graphitic systems the dominant energy contribution of x occurs when q and 



u are small so that we can approximate the bare response by its small q and u ex 



pansion 



2^ We 



xiq^u) ^ —{2h)~^q^{vlq^ + u^)"^/^ as calculated by DWR and Eqn. 3 of ref. 
can now write C = —k[1 + / (voq^]'^^"^ where n = j^— = 12.1 for graphene where 
Vq = 5.7 X lO^ms"-^. 

If we make changes of variables 6 = qD and sinh(?7) = ^ (so that C = —k/ cosh(r7)) then 
we can eliminate D from inside the integrals^. We thus obtain the energy expression 

t/vdw I' clA e'de cosh(r;)dry 

X [^a(^, v) + /t(cosh(r7) + XkY^] (8) 
= C^D-^ (9) 
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with 

J^x{e,r]) = — ^ ' (10) 

A/[cosh(^) cosh(r7) + XKsmh.{6)Y — cosh(?7)2 

and where the second term of equation ([S]) arises from letting D ^ oo in equation fllUI) 

Equation [8] is independent of D aside from the desired D^^ term so that C3 = D^Uvdw 
depends only on n. For the graphitic case where k = 12.1 we find 

C3 = 2.12 X 10^2 = O.SOeVA /atom. (11) 

By contrast the C4 coefficient predicted by Girifalco et al^ is C4 = 9.7949eVA^/atom which 
gives a potential approximately four times (0.079eV vs 0.022eV) as large as that of the 
inverse cubic power law at the experimental interlayer spacing Dq = 3.337A (equivalently 
this means that C^D-'^ > C^D-^ for D > ADq). 



III. NON-ASYMPTOTIC BEHAVIOUR 

While the C^D^^ power law will certainly be the dominant contributor to the dispersion 
potential for D ^ Dq, the intermediate-range (when D k, Dq) will include a number of other 
correlation effects. These include the C^^D^'^ potential from the atomic polarisibilities in the 
z direction in addition to a C^i2{D)D~^/'^ potential from the metallic electrons promoted 
from the vr^ orbitals due to layer overlap and hopping. As C5/2(-D) comes entirely from 
overlap of the vr^ orbitals it ought to be derivable from an analysis of the band-structure. It 
is expected to decay as an inverse exponential in D due to localisation of the vr^ orbitals. The 
C4 coefficient should be largely independent of D although some electrons will be promoted 
to metallic and graphitic response. 

With such a varied collection of correlation effects it seems unlikely that any simple ah 
initio method will adequately include the physics in the intermediate-range. Full RPA- 
ACFFDT calculations would be expected to provide a seamless potential through a wide- 
range of D however these are extremely difficult with current numerical approaches: for 
example, the van der Waals energetics of the semiconducting layered boron nitride sys- 
tem have been described succesfully using RPA energies,— but graphite gives convergence 
difficulties.— 
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IV. SEMIEMPIRICAL METHOD 



LDA calculations are expected to yield fairly accurate total energies for graphene when the 
interlayer spacing is compressed from its equilibrium value. Likewise the C^D~^ dispersion 
potential is expected to be accurate for layer spacing much greater than that of equilibrium. 
The intermediate range is more difficult to predict with neither method dealing sufficiently 
with the physics in that region. 

The method proposed in HN^ gives a fairly simple means (with minimal empirical con- 
tribution) of connecting the two regimes through the use of a fitting function. It is a 
semiempirical approach as the fitting function has its parameters chosen by matching exper- 
imental values for the lattice spacing and elastic constant C33. While this method predicts 
a reasonable value for the cohesion energy of graphite it, as with other methods, does not 
exhibit the correct behaviour in the tail due to the incorrect use of a C^D"'^ type dispersion 
law. In order to maintain consistency with this earlier work we re-examine the major results 
of their paper utilising the correct CsD"'^ dispersion law. To further maintain consistency 
we use the parametrisation of the LDA and GGA from the same paper. 



A. Semiempirical approach with pure C3D dispersion 

For our first new approach we adapt Equation 5 of HN to include the corrected form of 
the dispersion potential 

U{D) = [1 - f,{D)] Uj,MD) + UD)U,^^{D) (12) 

where [/vdw(-D) = C^D~^. Following HN, we use a Thomas-Fermi damping function 

/,(D) = [l + e-(^~^-)/^]-^ (13) 

where Dy/ and 5 are free parameters. The term involving A(^(4) is absent due to our C3 
coefficient being sourced from a bulk rather than a sum over pairwise potentials for multiple 
layers. U^^'^{D) is the parametrised LDA or GGA potential taken from equation 2 of HN. 

As in HN we attempted to determine 6 and Dy/ by ensuring that ^U{Do) = and 
^U{Do) = c-ss/ipDo) where 033 = 40.7GPa, Dq = 3.337A and p = 0.382A-2 take their 



experimental values (from ref. 26[ for C33 and ref. [271] for p and Dq). Using Uvdw = 
O.SOmeVA^D"^ we find that the HN fitting equations do not have a solution for the LDA 



or GGA. This lack of solution is not unexpected as the lack of other dispersion terms is 
expected to underestimate the dispersion for values of D ^ Dq. 

B. Semiempirical approach with mixed C^D^^ and C^D^^ dispersion 

While the CsD~^ term will certainly dominate over C^D"^ for D ^ Dq we know that it 
insufficiently models the physics for D ^ Dq, which we believe to be the cause of the fitting 
problems with the HN method for the semi-empirical method given in (A) above. The 
C4 = 9.795meVA^ coefficient used in HN is derived from a Cg = 16.34meVA^ coefficient 
calculated by Girifalco et al^ and constructed to ensure good Lennard- Jones modelling for 
a wide variety of graphitic systems. As such we propose to use its presumed accuracy for 
D Dq as a, correction to our C-^D^^ van der Waals function in order to better include the 
intermediate range physics. 

The simplest way to do this is to assume a correction to our function of the form C/^D''^. 
The C4 term is introduced firstly to cover the dispersion interaction due to the polarizability 
of the Hz and sp^ electrons in the z direction perpendicular to the graphene planes, and 
polarisability of the sp2 electrons parallel to the plane. These contributions to the dispersion 
interaction do not require long- wavelength collective electronic motions and therefor©^ are 
presumably describable by conventional asymptotics. These interaction are not included in 
our C^D^^ term, which is solely due to polarizability of the vr^ electrons along the graphene 
planes. The C4 term also has to account for the doped, metallic nature of the graphene 
planes near D ^ Dq due to overlapping of electron bands arising from the hopping of 
electrons from layer to layer. The corresponding attraction depends on the doping level, 
which decays exponentially with D. While this is not a law, it does decay faster than 
and hence is reasonably represented. 

Accordingly we now choose C4 so that the total van der Waals potential at the experi- 
mental lattice spacing remains the same in the two methods. This implies that 

C,D^^ = C^Df + C^D^^ (14) 

o 4 

which is true for C4 = 7.12meVA . This correction ensures we maintain similar D Dq 
behaviour while obtaining a correct power law for D ^ Dq. The van der Waals potential 
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now takes the form 



(15) 




FIG. 1: Potential energy versus lattice spacing (D). The solid line is the LDA corrected by the 
CsD~^ + C^^D^^ while the dashed line is the corrected GGA. The dash-dot and dotted lines are 
the pure LDA and GGA respectively. 

In order to ensure that equations [T2lfT3] correctly match the empirical data we must set 
5 = 0.221, Dw = 3.283 for the LDA and 6 = 0.340, Dw = 3.019 for the GGA when using 
the HN fitting function. Figure [T] shows the effect of this combined fit on both the LDA and 
GGA. 

In Figure [2] we show, for the LDA case, a more detailed comparison of three methods 
(the LDA, that of HN and the second method proposed here). It is quite clear from the 
graph that the method proposed here with the extra C*4D~^ correction closely matches that 
of HN for D ^ Do but maintains different asymptotes for D ^ Do- This suggests that the 
C^D^^ + C^^D^^ approximation, while a somewhat crude model of the true physics in the 
electron density overlap region, is able to maintain consistency with other methods. 

The most 'measurable' effect of the semiempirical approach is the interlayer cohesive 
energy min(f/(Z})). Table [T] summarises the results from HN with the addition of the new 
results calculated here. The new power law, used as the sole dispersion term does not give a 
valid cohesive energy due to the lack of a solution to the fitting function for both the LDA 
and GGA. The effect of the C^D^^ correction to the C^D^^ van der Waals functional on the 
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FIG. 2: Potential energy versus lattice spacing (D). The dashed line is the uncorrected LDA, the 
solid line is the LDA corrected by C^D^"^ + C4D~'^ while the dash-dots are the LDA corrected by 
C^D~^ as in HN. The inset shows the behaviour near Dq while the main graph shows the different 
asymptotics. 

cohesive energy is to give a very close cohesive energy to those predicted by HN, differing 
by only 1.7eV for the LDA and 2.3eV for the GGA or about 2 - 3%. 



LDA/GGA 


Expt" 


Expt'' 


Ucoh 26.5/2.3 


52.5 ±5 


or,+15 


L/G-vdW3 


L/G-vdW4 


L/G- 
vdW3+4 


Ucoh -/- 


60.7/57.4 


62.4/59.7 



TABLE L Cohesive energies per atom calculated by various approximations (in — meV). L/G- 
vdW3 is the LDA/GGA with a CsZ)"^ correction (for which meaningful results do not exist) while 
L/G-vdW4 has a dD-'^ correction (taken from HN using Cg = 16.34eVA^) and L/G-vdW3+4 



has the combined correction C-^D + C^D ^. The experimental results are taken from ref. 
for a and ref. j29] for h. 



28| 
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V. CONCLUSIONS AND FURTHER WORK 



In this paper we have investigated the asymptotic dispersion potential of 'stretched' 
graphite and found it to obey a C^D^^ type power law as opposed to the commonly employed 
C4D~'^. This places it in the same class of power law as bigraphene but in a different class to 
layered insulators (D~^ power laws) and layered metals (D"^^"^). This result has important 
implications for semiempirical (and otherwise) corrections to the LDA/GGA which have 
employed an incorrect power law. 

Furthermore we have employed the corrected power law in the simple semiempirical 
method of Hasegawa and Nishidate^ and found that, used as the sole dispersion term, it 
will not allow a valid solution to the HN fitting function. Reinclusion of a reduced C4D~^ 
term to include other physics from the non-asymptotic regime allows the method to be 
employed and gives similar results to HN for D Dq whilst ensuring the correct asymptotic 
behaviour is maintained. Its effect on the cohesive energy is fairly minimal with an absolute 
increase of the predicted cohesive energy of approximately 2 — 3%. 

While we believe that this power law (and the semi-empirical correction to it) should 
be accurate and useful for large layer spacings as in non-metallic intercalates and pillared 
systems we are not convinced that it will be as accurate in predicting the behaviour in the 
intermediate range of spacings without correction for other effects. Accurate RPA-ACFFDT 
calculations would provide a valuable benchmark for this and other methods. Until such 
time as these are available we hope that semi-emprical techniques like that discussed here 
should improve the accuracy of LDA and GGA calculations with widely spaced graphene 
layers. 
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